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Abstract 

In this paper we study the kinetic theory of many-particle astrophysical systems imposing 
axial symmetry and extending our previous analysis in Phys. Rev. D 83, 123007 (2011). 
Starting from a Newtonian model describing a collisionless self-gravitating gas, we develop 
a framework to include systematically the first general relativistic corrections to the matter 
distribution and gravitational potentials for general stationary systems. Then, we use our 
method to obtain particular solutions for the case of the Morgan & Morgan disks. The 
models obtained are fully analytical and correspond to the post-Newtonian generalizations of 
classical ones. We explore some properties of the models in order to estimate the importance 
of post-Newtonian corrections and we find that, contrary to the expectations, the main 
modifications appear far from the galaxy cores. As a by-product of this investigation we 
derive the corrected version of the tensor virial theorem. For stationary systems we recover 
the same result as in the Newtonian theory. However, for time dependent backgrounds we 
find that there is an extra piece that contributes to the variation of the inertia tensor. 
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1 Introduction 

The dynamics and evolution of collisionless stellar ensembles is a subject of great interest 
in astrophysics since they are the primary tool for comparisons of observations and the- 
ory in galactic dynamics pp. Such systems are generally composed by billions of stars, so 
it is neither practical nor worthwhile to follow the orbit of each particle in the ensemble. 
Most testable predictions depend only on the distribution function F(x, v,t) (DF), a quan- 
tity that determines the probability of finding a single star in a given phase-space volume 
<i 3 x<i 3 v around the position x and the velocity v. The dynamics of the DF follows from the 
appropriate kinetic equation and it in turn determines the statistical evolution of the system. 

In the framework of the general theory of relativity (GR) it is assumed that the DF satis- 
fies the general relativistic version of the Fokker-Planck equation [2j [3], H] or the collisionless 
Boltzmann equation (CBE) E]. The first one is devoted to systems in which local gravi- 
tational encounters dominate in their evolution whereas the latter is useful to study systems 
sufficiently smooth, so that they may be considered to be collisionless pQ. One can actually 
consider systems in which a number of particle species can collide and produce different 
species. This is how the formation of the light elements in the big bang nucleosynthesis is 
calculated (see [7] for a review). However, in systems such as galaxies and galaxy clusters, 
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physical encounters between the stars are very rare, and the effect of gravitational collisions 
can be neglected for times far longer than the age of the universe. Thus, for these systems, 
the CBE provides a very good approximation. 

For a typical galaxy, the relaxation time, t re \ ax , is arbitrarily large in comparison with 
the crossing time, t CT oss- This means that they can be approximated as a continuum rather 
than concentrated into nearly pointlike stars. Now, it is commonly assumed that the main 
contribution of the mass in a galaxy is concentrated in an axisymmetric flat distribution 
PQ. For this reason the obtention of idealized thin disk models has been a problem of 
great astrophysical relevance. In this case, the most straightforward way to construct a 
self-consistent model is by means of finding the DF for a system with a known gravitational 
interactions and matter distribution. Since the mass density of the system is defined by the 
integration of the DF over the velocity, the problem of finding a DF is that of solving an 
integral equation (see [51 13 EDI Ell EE21 EES] and the references therein). At present we have 
at disposal a variety of self- consistent galaxy models: [HI 13 HH HH HH UHl [201 EU [22l E31 

12a lis]. 

Now, even though for most systems under consideration Newtonian gravity is believed to 
be dominant, general relativistic corrections might play an important role in their evolution. 
As a matter of fact, in recent years it has been an increasing interest in the incorporation 
of GR in the description of these systems, and to date we have a variety of fully relativistic 
galaxy models: ^SII2ZII2HII221EDIEIII321E31I33, among others. Perhaps the principal 
reason of including GR corrections in galactic dynamics, is the hypothesis that it is possible 
to overcome the problem of the rotation curves predicted by the Newtonian theory. While, 
some authors argue that by using GR the inclusion of a dark matter halo is unnecessary at 
galactic scales (see for instance [2S1 ESI EZ1 ESI ES]), several publications have pointed out 
that this is not entirely true [101 HU H21 [331 E5J . In particular, the authors of [35] presented 
a model in which the percentage of dark matter needed to explain flat rotation curves turns 
out to be ~30% less than the required by the Newtonian theoryfj]. It is important to point 
out that currently there are alternative approaches to GR which address the problem of 
rotation curves in spiral galaxies, as for example the so-called f(R) gravity (see [47J, for 
references) . 

Despite the fact that the relativistic contributions do not solve completely the problem 
of rotation curves in galaxies, it seems that they do introduce significant corrections. Thus, 
in order to estimate the effects on the various observables we are interested in, it would 
be nice to have a framework to include systematically general relativistic corrections to 
a given Newtonian model. The post-Newtonian approximation is perfectly suited for this 
purpose. The appropriate scheme that describes the effects of the first corrections beyond 
the Newtonian theory, was first formulated in [381 1391 [50] (see [51] for a textbook analysis) 
and it is known as the first post-Newtonian (1PN) approximation. This approach holds 
if the particles in the system are moving nonrelativistically (v <C c) (as in the case of a 
star moving around a typical galaxy) and gives the corrections up to order v 2 /c 2 , where v 
is a typical velocity in the system and c is the speed of light. Currently, higher order PN 

1 Strictly speaking, the aim of introducing a dark matter halo is not necessarily to reproduce a flat rotation 
curve but generate a curve that rises and then declines at particular paces [46] . 
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approximations have been developed because of the increasing interest around kinematics 
and associated emission of gravitational waves by binary pulsars, neutron stars and black 
holes, with promising candidates for detectors such as LIGO, VIRGO and GEO600 (see 
[52] 153] for references) . 

Based in the above considerations, we recently started a general study of self-gravitating 
gases in the collisionless regime and, as a first step, we derived a version of the CBE that 
accounts for the first general relativistic corrections [51] . With this tool in hand, we obtained 
the 1PN version of the Eddington's polytropes, starting from an ergodic DF proportional to 
E n . The purpose of this paper is twofold. First, to implement a similar procedure in the 
axially symmetric case in order to setup our general framework. And second, to obtain a 
new set of self-consistent models starting from a Newtonian "seed" and study the impact of 
relativistic corrections on the various observables. 

The rest of the paper is organized as follows. In section [2] we present a brief overview 
about the basics of the 1PN approximation, revisiting the field equations, as well as the 
kinetic theory for arbitrary self-gravitating systems. In section [3] we show the fundamental 
equations defining self-consistent models with post-Newtonian corrections. We start dealing 
with arbitrary systems but then we focus on discoidal configurations with axial symmetry, 
in order to prepare the ground to construct 1PN galaxy models in section HI Finally, we 
summarize the principal results in section |5j 

2 General Framework 
2.1 The 1PN approximation 

The post-Newtonian approximation has been reviewed carefully in a number of references 
(see for example [ST]). However, we will include here the basic definitions and relations for 
completeness. 

First off, note that in Newtonian mechanics the typical kinetic energy is roughly of the 
same order of magnitude as the typical potential energy, and thus 

v 2 ~ 4>, 

where v is the mean velocity in the system. The idea is then to express all physical quantities 
in terms of a series expansion of the small parameter i)/c C 1 and keep the leading order 
beyond the Newtonian theory. The first quantity to consider is the spacetime itself: any ma- 
nifold can be considered to be locally flat so, for particles that are moving nonrelativistically, 
we proceed to express the metric tensor as 

2 4 

900 = -1+ 9oo + fl'oo H , 

2 4 

9ij = 9ij + 9ij H , (1) 

13 5 

901 = 90i + 9(H + 90i + " " " j 
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where the symbol g^ u denotes the term in g^ v of order (v/c) n . Odd powers of v/c appear in 
goi because these components must change sign under time-reversal transformation t — > —t. 

It is natural to assume a similar expansion for the components of the energy momentum 
tensor. From their interpretation as the energy density, momentum flux, and energy flux, 
we expect that 



o 

X 00 = T 00 + T 00 + • • • 

2 4 

T ij = T ij + T ij + . . . ^ ( 2 ) 

1 3 

rpOi rpQi rpOi 

These expansions lead to a consistent solution of Einstein field equations. 

Working in harmonic coordinates (i.e., coordinates such that g^T^ = 0) and to our 
order of approximation, the various components of the metric tensor can be expressed in 
terms of the Newtonian potential and post-Newtonian potentials ip and & as 

#00 = -20/c 2 , 
Joo= -2(0 2 + ^)/c 4 , 

2 9ij = -20%/c 2 , (3) 



9oi = 0, 

3 



goi = £i/c . 

Thus, the Einstein equations reduce to 



o 

V 2 = AuG T 00 , (4) 



2 2 Q2, 

72„/, _ A„n^i( T-00 



V 2 ^ = 47rGc 2 ( T uu + T M ) + — I , (5) 



1 

V 2 & = 167rGc T°\ (6) 
along with the coordinate condition 

4^ + V-£ = 0. (7) 

One can also consider the motion of test particles in a given background. For general 
potentials 0, ip, and ^ one finds that the free falling particle obeys the equation 



V (20 2 + f)+^-vx(Vx^)- 3v^ - 4v(v • V0) + w 2 V0 



(8) 



which partially resembles the mathematical structure of the Lorentz force experienced by 
a charged particle, with velocity v, in the presence of an electromagnetic field. Such law 
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of motion will determine, for instance, the rotation curve corresponding to a given galactic 
model (see Appendix IA1 for details). 

It is instructive to point out that the equations of motion (|S]) can be derived from the 
Lagrangian (51] 



2= V - 
2 



■7 T + 



— \~ 1p 



(9) 



For stationary spacetimes, the potentials are independent of time and the associated Hamil- 
tonian % = ^2 i Xid£ / /dVi — £, is a conserved quantity that can be interpreted as the 1PN 
generalization of the classical energy: 



E 



3v A 3v 2 



(10) 



Note that this expression is independent of the vector field 

If the source of gravitation is endowed with axial symmetry, the z-component of the 
angular momentum is an additional integral of motion. In cylindrical coordinates (R, if, z) 
we obtain that, for (^-independent potentials, the quantity 



RV<n + 



Rv u 



V 

y 



30 + RL 



can be interpreted as the 1PN generalization of the azimuthal angular momentum. In this 
case £ plays a role, through its rotational component. 



2.2 1PN statistical mechanics 

From a statistical point of view, the state of the system can be determined by its DF, 
F(x, v,t), depending on the spatial coordinates, velocity and time. Now, as mentioned in 
the Introduction, for applications in galactic dynamics it is commonly assumed that the 
encounters between particles are negligible and hence, the evolution of the stellar system 
must obey the so-called collisionless Boltzmann equation. In 1PN approximation, such 
relation can be written as [54J 

dF { dF d(j)dF 1 / V ,\ (9F { dF\ 



+1 L V ^ - (— + 3^1 ^ -v j - ?k) +3v^ -^-^1^=0 
c 2 dxi \ 2 J dx l \dxi dx l J dt dx l dt dv l 

(12) 

For situations in which encounters play a dominant role, the right-hand side of the above 
expression can be replaced by a collisional term of the Fokker-Planck type [55] . 

If the self-gravitating system is in a stationary state, dF/dt = and the DF can be 
expressed as a function of the energy (a first integral of motion). Moreover, if the system 
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is endowed with additional symmetries we can expect that, as in Newtonian theory, the 
stationary solutions of ( I12p are also functions of the remaining integrals of motion. In other 
words, there is a 1PN version of the Jeans theorem, which simply reflects the fact that 
equation (fl~2l can be rewritten as dF/ dt = [54J . 

The DF describes completely the sate of the system. We can extract from it all the 
relevant physical quantities, like the mass density, the mean velocity, the velocity-dispersion 
tensor, etc. For instance, by taking the moments of equation f[T2"j) . we could in principle 
obtain some useful results regarding the hydrodynamics of self-gravitating systems at 1PN 
order. This was first derived by Chandrasekhar [56J using a different approach so we will 
refrain from writing out these results here, since they are not particularly illuminating. 

Another important fact that can be derived with the help of the DF, is the 1PN virial 
theorem (see appendix|B]for details). This theorem shows that there exists a linear relation 
between the variation of the moment of inertia tensor P k , the kinetic energy tensor K lk and 
the potential energy tensor W , according to: 



This 1PN virial equation has a similar form as the Newtonian one, except for the last term 
in the right-hand side. The temporal variation of <fi contributes to the variation of I . In 
particular, for stationary states we have 



as in the Newtonian case. This is, the total gravitational potential energy (a negative 
quantity), is two times the total kinetic energy. 

3 The 1PN Self-gravitation Equations 

The integration of the DF over the velocity space, leads to the energy-momentum tensor. It 
in turn determines the gravitational field via the Einstein equations, so it should be possible 
to derive a relation between F and the potentials <ft, ip, and Such relation is represented 
by a set of couple equations which we shall call the 1PN self- gravitation equations. At first 
we deal with the general case of arbitrary self-gravitating systems and then we shall consider 
the special case of razor thin disks. 

3.1 Fundamental equations for arbitrary systems 

We begin by considering the following expression for the energy-momentum tensor, valid for 
any self-gravitating system: 




o 



(13) 



W 



-2K. 




(14) 
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Here U v is the four- velocity, which is related to the three- velocity by U l = U°v l /c (Greek 
indices run from to 3 and Latin indices from 1 to 3). From the relation g ilv U il U v = —c 2 
we can compute the quantities we need to derive the 1PN approximation of the energy- 
momentum tensor. After some calculations, we obtain 

- = l+ 2*-^ (15) 



and 

d 3 TI =114- 

2c 2 c 2 

Now, taking into account that \f—g = 1 — 20/c 2 + we find that 



T 00 = / -fFd 3 v, 



T j = i / 7 vVFd 3 t;, (17) 
T 0i = - c I - r'Fd'r. 



where 



3v 



2 



T = 1 + ~2 2 

c c 



is the measure of the integration over velocities. 

According to the above relations, we expect that the DF can be expanded in power series 
of v/c as 

F = F+ ! + •••, (19) 

o 

where F is the Newtonian contribution to the DF (which is itself a solution of the classical 

2 

CBE) and F is the first post-Newtonian correction. Plugging (fl9l) into ( FlTl) leads to the 
different components of the energy-momentum tensor at the orders required by the 1PN 
approximation: 



o 



T 00 = / Fd 3 v, (20) 



2 



T 00 = -| / (v 2 - 20) F d 3 v + / J 1 rf 3 w, (21) 
2 1 /" o 

T l > = - v L v 3 F d 3 v, (22) 



c 2 



i 



T 0i = - [ v l F d s v, (23) 
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along with T* J = 0, as expected. With these results we are ready to write the 1PN self- 
gravitation equations for general stationary systems: 

V 2 = 4ttG J Fd 3 v, (24) 
VV = 8ttG / (2v 2 - 30) F d 3 v 

.2 / I .!■'> 



+ 4nGc 2 J Fd 6 v, (25) 

V 2 & = 16ttG y F (26) 

To summarize, we can say that a stellar system characterized by an equilibrium DF satisfying 
(]I2"]) is described by a matter distribution given by (12"U]) - (12"3"]) . and gravitational interactions 
determined by the field equations (jH)-(jSD- In order to provide a self-consistent description, the 
relations (|24j) - (|26p must be satisfied. All of these equations are written as power expansions 
in the small parameter v/c and, in consequence, we can clearly distinguish between the 
Newtonian contribution and the post-Newtonian corrections. 

3.2 The case of stationary razor thin disks 

The so-called razor thin disks are of special interest in modeling a number of axisymmetric 
galaxies. In this case, the DF depends on velocities and positions in the form 

F = f(R,v R ,v ip )6(z)6(v g ), (27) 

where 5 is the Dirac delta function and / is a reduced phase-space density describing the 
stellar population placed on the equatorial plane z = (for a finite thin disk of radius a, / 
must vanish for R> a). Equations (|24|) and (125]) can be written as 

V 2 = 4rrG5(z) [ f d 2 v, (28) 



VV = ^G8{z) 



(Av 2 - 60) f d 2 v + c 2 [ f d 2 v 



(29) 

where d 2 v = dvftdv^. Equation (126]) requires a little more attention. First start with equation 
(JSJ) (which is equivalent to (T26]) ). The general solution that vanishes at infinity can be written 
as 0] 

i 

6(x) = -4Gc / T °* WV , (30) 
J l x ~ x l 

which, in our case, reduces to 

-4G f [ - - d3x ' d3v AG f 5 ^ z ') d ^ x ' d2v 
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where we have computed the integral with respect to v z (we use the notation v 1 = v x , v 2 = v y 
and v 3 = v z ). This expression can be massaged into a useful form by taking into account 
two facts: (i) the relation between the cartesian components v x , v y and the cylindrical 
components vr, v v , i.e. v x = VRCOsf — v^simp and v y = VRsimp + v^cosf, (ii) since we 
are dealing with stationary axisymmetric systems, the DF is an even function of vr [1] and 
in consequence J vr/cIvr = . Thus, we can write 

v x f 5(z')d 3 x'd 2 v f f v v f 5(z')d 3 x'd 2 v 

: : = — Sm<p 





and 

v y f 5{z')d 3 x'd 2 v _ cos(fi f fv v f 5{z')d 3 x'd 2 v 





x — X' 



Now, by introducing the relations £ x = £r cos if — ^ sin ip and £ y = £r sin ip +£ v cos (p in ([30 
we obtain 

o 

v v f S(z')d 3 x'd 2 v 



£r cos (p — sin if — sin ip / / AG 

J J I* - * | 

and 

t ■ ,c f f A ^ v/f5(z')d 3 x'd 2 t 
(Hsm^ + L cos f = — cos f 4G — - — ; 



Since we assume that £ is ^-independent, each of these expressions leads us to the conclusion 
that 

£ij = 0, for finite distributions, (31) 
and that ^ is solution of the following equation: 

V% = 16ttG5(z) fv v f d 2 v. (32) 

The equation for the component £ z can be obtained easily by replacing ( |2"7|) in ([26]) . and the 
result is the Laplace equation, V 2 £ 2 = 0. Its solution can be determined through condition 
(Jjj. A straightforward calculation leads to 

UR)=UHR/Ro), (33) 

where £ zo and R Q are constants of integration. In the case of distributions with finite extent, 
we demand as a boundary condition that lim^oo £ z = 0. In consequence, we have to choose 
£ zo = 0, and hence 

£ z = 0, for finite thin disks. (34) 

On the other hand, we expect that / obeys the collisionless Boltzmann equation in the 
three-dimensional phase-space (R, vr, v v ). In fact, by introducing ( 1271) in ( fl2|) and performing 
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an integration on z and v z , it follows that the distribution / obeys the relation 




c 2 c 2 dR 



<f) 4Rd(j) 



c 2 v v dR J R dv, 



R_diA vrv<p df 




(35) 



+ 




+ 



3v 2 — 5v i 



2c 2 



c 2 7 <9# c 2 aR c 2 9/2 J dv R 




where v 2 = v 2 R + t> 2 and d<fi/dR, dip/dR and d^/dR are evaluated at z = 0. The above 
relation is the 1PN version of the Boltzmann equation for an axisymmetric two-dimensional 
shell, located at the equatorial plane, and in a stationary state. Of course, f(R, Vr,v v ) plays 
the role of the reduced DF describing the diskoidal shell. 

It is straightforward to show that E and L z , given by equations flTUJ) and ( ITT]) , are 
solutions of 035]). This means that, for axially symmetric systems, any / depending on E 
and L z is solution of the CBE; conversely, any solution of the CBE can always be expressed as 
a function of E and L z . Thus, any two-integral DF, f(E, L z ), provides a complete statistical 
description for the (two-degree-of-freedom) stellar system. This fact will be very useful for 
the formulation of post-Newtonian models in the next section. 

4 Analytical Models for Axisymmetric Galaxies 

The purpose of this section is to show how to implement the formalism developed above in 
order to obtain axially symmetric galaxy models. For the applications we want to consider 
here we have to take into account further considerations. First of all, recall that in [51] we 
proved that Jeans theorem remains valid at 1PN order. This means that any equilibrium 
solution of the CBE depends only on the integrals of motion of the system, and that any 
function of the integrals yields an equilibrium solution of the CBE. Thus, for stationary 
systems with axial symmetry, we can restrict ourselves to DFs depending on the energy ffTU]) 
and the angular momentum ( ITTj) . which are themselves integrals of (Il2p. 

The next step would be to implement the previous restrictions starting from a given 
Newtonian potential-density pair with a known DF, as was done in [53] for the spherically 
symmetric case. As a result, one expects two coupled self- gravitation equations, providing a 
method to determine, from a Newtonian model, its associated post-Newtonian corrections. 
In practice, the present formalism leads to a two coupled ordinary differential equations 
in the spherically symmetric case. In the axially symmetric situation however, one ends 
up with two coupled elliptic partial differential equations (for general volumetric matter 
distributions). 

In this case such equations are much more involved than the ones corresponding to the 
spherically symmetric case, but the configurations we shall deal here permit us to introduce 
some additional assumptions to simplify the problem. In the next section, we shall show that 
a dramatic simplification can be achieved by the consideration of thin discoidal distributions 
in spheroidal oblate coordinates: instead of getting differential equations, in this case the 
post-Newtonian corrections can be obtained from simple algebraic equations. 
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We then present a particular application where the resulting equations can be solved 
analytically, which means that it is possible to obtain 1PN exact solutions. The importance 
of these solutions will be evaluated by a comparison between density profiles and rotation 
curves described by Newtonian theory and the ones predicted by the 1PN approximation. 
Although focus on the particular models introduced in [19] (revisited by [21]), our framework 
can be applied to a wider variety of models. In general, the method can be used for situations 
in which the potentials are separable functions of the spheroidal oblate coordinates. 



4.1 Hunter's method in the 1PN approximation 

One can find in the literature a number of self-consistent stellar models representing razor 
thin disks; here we will deal with models belonging to the family of Morgan & Morgan disks 
[T9"[ [21]. In the Newtonian formulation, they can be obtained by a formalism developed by 
Hunter [57] . Such procedure (known as the Hunter's method) provides the surface density 
of the disks, the gravitational potential, and the circular velocity as a series of elementary 
functions, by superposing solutions of the Laplace equation in oblate spheroidal coordinates. 
Hunter's method can also be implemented in the context of the 1PN approximation as 
follows. 

To begin with, note that in vacuum, the field equations flU)-© reduce to three Laplace 
equations for 0, ip, and £ v (remember that £ 2 = for distributions of finite extent). Without 
loss of generality, we assume that the disk is on the equatorial plane, so we have to impose 
that the gravitational potentials have symmetry of reflection with respect to the plane z = 0, 
i.e., <j>(R,z) = 4>(R,—z), if)(R,z) = i/)(R,—z), and £ v (R,z) = £ v (R,—z). Then, it follows 
that 

*(*,-,) = -*(*,), (36) 
»*(*,-,)=-§*(*,„, (37) 

ir«-*> = -il«*>. (38) 

in agreement with the attractive character of gravitation. We also assume that d(f)/dz, 
dip/dz, and d^/dz do not vanish in the disk's zone, in order to have the corresponding thin 
distribution of energy-momentum. Such distribution, restricted to a region < R < a in the 
plane z = (from here on, a will denote the disk radius), will be described by a "shell-like" 
energy-momentum tensor. If we define 

T 00 = E(R)6(z), (39) 



2 2 
~>00 i rpii 



1 



T w +T n = -a(R)5(z), (40) 



c 2 



T°* = -A(R)5(z), (41) 



c 



11 



for < R < a, it follows from Gauss's Law that 



-^) = ^(-) . (12) 

=0+ 



2ttG V dz 



— (%) (in 



=0+ 



Note that E represents the surface mass density of the Newtonian theory (i.e. without 
relativistic corrections), while A plays the role of the surface density of ^-momentum. On 
the other hand, a is associated both to the pressure and the relativistic corrections to the 
mass surface density. 

The above relations mean that, in order to have a distribution of matter as the described 
by (}4"2"]) - (j4"4"]) . we have to demand that 

^(R,0 + )^0, R<a, (45) 
oz 

^( J R,0 + ) = 0, R>a } (46) 
oz 

with the same requirement for if> and At this point it is convenient to introduce oblate 
spheroidal coordinates, a system that adapts in a natural way to the geometry of the problem. 
They are related to the cylindrical ones through 

R = a> /(l + C 2 )(W), (47) 
z = a(ri, (48) 

where < £ < oo and — 1 < rj < 1. Note that (i) the disk itself has coordinates ( = 0, 
rf = 1 — R 2 /a 2 ; (ii) conditions (j4"5]) - (j4"6"]l become 



dr] 



= H( V ), (49) 

C=o 

= 0, (50) 

J7=0 



where H is an even function of 77. The general solution of Laplace's equation satisfying the 
above conditions can be written as 

00 

<P(C,V) = -^A 2n g 2 „(C)P 2ri ,(r/), (51) 

n=0 

where A 2n are arbitrary constants, P^nijf) an d l2n(C) — i 2n+1 Q2n(i() are t ne usual Legendre 
polynomials and the Legendre functions of second kind, respectively. The post-Newtonian 
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potentials ip and £ v have the same form, 

oo 

1>{C,V) = -J2 B ^^P2n(v), (52) 

n=0 

oo 

UCV) = -J2 C ^(0P2n(v), (53) 

n=0 

but here we have denoted the expansion constants as B 2n and C 2n . We can derive explicit 
formulas for S, a, and A in oblate spheroidal coordinates, by introducing (}5"Tl) - (l5"3l in (|42|) - 



S = A ^ 2n + 1 )92n+l(0)P 2 n(^), (54) 

^ ^* n=0 
- x 

a= o ^ y]^2n(2n+l)g 2 » +1 (0)P 2 »(?7,), (55) 



n=0 



4 



— VC 2n (2n + l)g 2n+1 (0)F 2n ( )?t ), (56) 



n=0 



where 77* represents the value of coordinate 77 inside the disk: 



(57) 

Now that we have stated the fundamental structure of models with 1PN corrections, the 
next step is to demand that the models obtained are self-consistent, i.e., that they have an 
analytical equilibrium DF that is related consistently to the surface mass distribution. In 
other words, we have to formulate the corresponding 1PN self-gravitating equations: 

±A 2n ^hl = r ifv< (58 ) 

n=0 V * J 

B 2n ^^ = [ (Av 2 - 60) fcfv + c 2 [f d\ 

(59) 

jrCin^^- = [v/fd\ (60) 

n=0 ^ J 

where, for the sake of simplicity, we have defined 

~ (2n + l)g 2n+1 (0) 

A 2n = — A 2n , (61) 
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the same for B 2n , and 

~ _ (27i + l)g 2n+1 (0) ^ 

^2n — ~ -p, ^2n- 

8naG 

Note that, part of the right-hand side of equations (l58l) - (!60j) is characterized by three fun- 
damental quantities in Newtonian dynamics of self-gravitating systems: (i) the mass surface 

o o 
density, £ = J f d 2 v; (ii) the mean square velocity, (v 2 ) = E^ 1 J v 2 f d 2 v, and (iii) the 

o 

mean circular velocity, (v v ) = E~ x J v v f d 2 v. In particular, equation ( 158]) is the Newtonian 
self-gravitation equation. 

There is a variety of cases that can be addressed by the formalism presented above. They 
correspond to situations in which the gravitational fields (i) have symmetry of reflection 
with respect to the equatorial plane, (ii) are separable functions of the spheroidal oblate 
coordinates and (iii) correspond to axially symmetric distributions with finite extent. Note 
that any function with these features can be expressed in the form fl51R . 



4.2 The Morgan & Morgan disks with 1PN corrections 

In Newtonian gravity, the Generalized Kalnajs Disks [23] are finite thin distributions of 
matter with surface mass density 



(2m + l)M / _E?\ m - l/2 
2yra 2 V a 2 



E m (R) = 1 n _J ( 1 - — ) , m e M (63) 



and gravitational potential <p m , given by (I5TT) . with the following expansion constants 

(m) _ MG v ^2- 2m - 1 (4n + l)(2m+l)! 
^Sn — — — 77; TTT7 7T7v ; 77v\ • I 04 J 



These solutions were first obtained by Morgan & Morgan (MM) [19J, by solving a boundary 
value problem in the context of GR. Even so, we refer to the above solutions as MM disks 
for historical reasons, but keeping in mind that we are dealing with Newtonian gravity. 

All the members of this family have interesting features: a monotonically decreasing mass 
density and a Keplerian rotation curve. The case m = 1, corresponding to the well-known 
Kalnajs disk [20] (see also [TJ), is an exception because it describes a self-gravitating disk 
which rotates as a rigid body. 

The superposition of members belonging to the MM family, has been used to obtain 
new models with more realistic properties. For example, in reference [25], the authors con- 
structed a family of models with approximately flat rotation curves considering a particular 
combination of MM disks. Another example is the family of models introduced by Letelier 
in [58], representing flat rings (see also [59] for astrophysical applications). Additionally, in 
reference [60] it was shown that it is possible to construct models with realistic rotational 
curves obeying simple polynomial expressions. In particular, the authors constructed models 
for a number of galaxies of the Ursa Major cluster, and as an application they estimated 
their corresponding mass distributions. This fact suggest that there exist superpositions 
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of MM members leading to galactic models in agreement with the so-called maximum disk 
hypothesis [JJ. It would be interesting to have in hand the 1PN version of the MM family, 
along with all of its features. 

Here we illustrate how to obtain the post-Newtonian version of the MM disks. We focus 
in the m = 2 model, which is characterized by a density 



E 



5M 
2vra 2 



(65) 



but in principle, the same procedure can be applied to any of the remaining members (or 
linear combinations). From here on, we will drop the subindex of S m , for simplicity. 
It can be shown that the above surface density can be obtained from the DF [BT] 



where 



*=4= 



10M 



k{QL z - E - 5 n 2 a 2 /A)~ l/ \ 



y/3 [7r u a 5 G\ 



Q 



157rGM 
32a 3 



(66) 



(67) 



Note that the DF defined by ( J66l) is function of the Jacobi's integral, Ej = E — f2L z , which 
can be interpreted as the energy measured from a frame that rotates with constant angular 
speed Q (see Appendix ICl) . 

In order to obtain a 1PN version of the model, we start from the DF given by f )66|) but 
using the post-Newtonian expressions for E and L z , i.e., equations ( ITU]) and ( ITT]) . Thus, at 
1PN order we can write 

2 
/=/ + /, 



with 



where 



/= k,J 1/4 and /= -7^2^, 



Jn 



J-2 



— c 



v 2 /2 + f2Rv v - 5i?V/4, 

+ 1/J 7T- + — ORv u 



-5/4 



(68) 
(69) 

(70) 
(71) 



2 2 

These relations determine the 1PN self-gravitation equations through ( I58|) -(l60l) and the in 
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tegrals 



where 



fd 2 v 



v 2 f d 2 v 



f d 2 v 



fd 2 v 



2325 
1024 : 



5M 



2na? 

75GM 2 
448a 3 

5 /15GM 3 



c 2 7/* 



27ra 5 



(7_ 7 ^ + 6 ^ ^3 



3=0 



160M 

37r 3 a 3 G 



75 
512' 



1125 
128 : 



1275 
512 ' 



32^ 
32aGn 2 



38475 
" 7168 ' 



(72) 
(73) 
(74) 
(75) 



(76) 



By introducing (}72"]) - (j75]) in (j5"81) - (j6"01) . we obtain a system of linear equations for the constants 
B 2n and C 2n in terms of A 2n (which are known a priori), that can be solved analytically. 
After some computations we find that 



a 



( GM\ 



3/2 



2/i 



(45 - 48n - 92n 2 + 16n 3 + 16n 4 ) {An + 1) 5tt 5/2 ^15 



\ « j ra-n)ra-n)r(l + n)r(4 + n)(2n + l)q 2n+1 (0) 256^ 



(77) 



B 



2n 



GM 



.2 I V2 

l2n ~ Xlm=0 ^2n,2m 



a J (64/3)g 2n (0) + vr(2n + l)g 2n+1 (0) 



(7E 



where J 2n and B 2n ^ 2m are constants defined in the Appendix O, equations (j!24p and (j!26p . 
respectively 

With the constants C 2n and B 2n at hand, the post-Newtonian fields are completely de- 
termined along with the remaining physical quantities. In galactic dynamics, the circular 
velocity and the mass profile are two important measurable quantities used to verify a par- 
ticular model. In our case, the first one is given by the relation f )86|) whereas the second one 

2 

is given by (165]) plus the 1PN corrections coming from T 00 . The corrected surface density 
can be written explicitly as 



2 
E + E-- 



5M 
2na 2 



V* 



15?rA . 2 4 v 3 A {2n)\\B 2n 



224 



P2n{V*) 



(79) 



where A is a dimensionless parameter defined by 

GM 



A 



SO) 
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Figure 1: Rotation curves for the 1PN version of the second MM member for different values 
of the parameter A: (a) A = HT 3 , (b) A = 5 x 1(T 3 , (c) A = 2 x 1(T 2 , and (d) A = 7 x lO" 2 . 
For A ~ 10~ 2 or more, the differences with the Newtonian model (dashed line) become 
significant. In this case, the rotation curve tends to flatten as A grows. 

and B 2n are the dimensionless constants 

B2n = (81) 

The parameter A, which is also present in the expression for the circular velocity, is a 
measure of how large are the 1PN corrections. For example, in a galaxy with 10 12 solar 
masses and a radius of 10 Kpc, we have A ~ 5 x 10~ 6 . Here, we consider situations where 
the relativistic corrections are larger and they can be visualized in the behavior of rotation 
curves and mass profile. In Figure (TJ we plot the circular velocity when A ~ 10~ 3 and 
A ~ 10 -2 . The 1PN corrections become important in the latter case (last two figures), in 
particular for values of v v near to the disk edge. This is somewhat surprising since one would 
expect major corrections near the galaxy core, where the mass concentration is maximum. 
A similar phenomenon occurs with the mass profiles (see Figure [2]). Their differences with 
the Newtonian profile become significant for A ~ 10 -2 and the magnitude of 1PN corrections 
increase with the radius. 



17 




R/a R/a 



Figure 2: Mass profiles for the 1PN version of the second MM member for the same values 
of the figure [lj As expected, the differences with the Newtonian model (dashed line) become 
significant for A ~ 10~ 2 or more. In the case (d), we have a profile in which the maximum 
is not at R = 0. 



5 Concluding Remarks 

We continued our study of the post-Newtonian kinetic theory of collisionless self-gravitating 
gases, extending our previous results [54J to systems with axial symmetry. Before dealing 
with the applications, we developed the 1PN version of the tensor virial theorem, applicable 
for arbitrary self-gravitating systems. We found that the 1PN virial theorem differs from the 
Newtonian one by the fact that the temporal variation of <fi contributes to the variation of the 
inertia tensor. However, for stationary systems we recover the same result as in Newtonian 
theory: the absolute value of the total gravitational potential energy is two times the kinetic 
energy. The problem about the validity of the virial theorem for collisional systems (which 
remains valid in the Newtonian theory) is an open question that we leave for future studies. 

The applications we considered here were devoted specially to address the modeling 
of systems in galactic dynamics, thus we focused in equilibrium situations with rotational 
symmetry. We developed a formalism to construct the 1PN version of Newtonian self- 
consistent models characterized by a stationary DF depending on the Jacobi's integral. The 
general case is mathematically challenging as it involves a pair of coupled elliptic partial 
differential equations. However, we restricted our attention to the study of thin models with 
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finite extent, the case in which the system reduces to a set of algebraic equations. As an 
illustrative example of the formalism developed here, we obtained the 1PN version of the 
second MM disk. The contributions of relativistic corrections found are not significant for 
A ~ 10~ 3 or less, but they become prominent for larger values of A (provided that the 1PN 
approximation is still valid). We noted that the corrections grow with the radial distance, 
which is surprising since one would expect major corrections near the disk center, where the 
mass concentration is maximum. 

The general formalism starts by implementing the Hunter's method in the 1PN scheme. 
Then, a stationary DF (solution of the 1PN CBE) is selected and introduced in the self- 
gravitation equations. Finally, by solving the resulting algebraic equalities, we find the 
corresponding post-Newtonian fields, if) and The new lPN-corrected models obtained 
are also self-consistent, since this feature is an inherent requirement of the formalism. The 
method is applicable to self-gravitating thin disks with finite extent provided that: (i) the 
system has symmetry of reflection with respect to the z = plane; (ii) the Newtonian po- 
tential is separable in spheroidal oblate coordinates; (iii) The Newtonian DF is only function 
of E and L z . 

It is worthwhile to point out the interesting work by Schenk, Shapiro and Teukolsky in 
1999 [52] ■ In this paper the authors focused on the Kalnajs disk (the m — 1 member of the 
MM family) and solved numerically the Einstein equations coupled to the relativistic CBE, 
obtaining the first self-consistent model of a rotating relativistic disk. It would be interesting 
to obtain the 1PN version of the same model by using the formalism developed in this paper 
in order to establish a comparison of the two approaches and estimate the effects of the 
higher order relativistic corrections. 

As a final remark, it is important to mention that the problem about the stability of 
these models is an important subject that will be addressed in future works. This issue 
brings several points to be investigated. One fundamental question is to establish whether 
the discoidal structures considered here are stable by themselves or if they need an additional 
component (a spheroidal halo, for example) in order to be stable. Another question is to 
determinate if the marginal stability of Newtonian disks improves or not with the introduc- 
tion of relativistic corrections. At the present we are performing a preliminary study about 
the stability of test particles around the solutions obtained, and their response to radial and 
vertical perturbations in order to compare the results with the Newtonian theory (see for 
example [211 E3]) and to provide a first test of stability. Then we expect to implement a 
more conclusive analysis regarding the Jeans-type instabilities by perturbing the DF's by 
themselves. 
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A Rotation Curves in the 1PN Approximation 

The circular velocity is a widely used observable in galactic dynamics, so it is important to 
have in hand a formula for the rotation curves in the 1PN approximation. In order to do 
this, we have to study the circular motion of test particles in the equatorial plane. At first, 
note that according to (jSJ), a star moving around an axisymmetric stationary thin disk with 
finite extension, obeys the following equations of motion: 

R = R ^_^±( 1 + ^U (j) _ 3R 2 + R 2 (p 2 + ^ 



dR 



+ n iRi t-m +R m)> (82) 



dz \ 6 

Rip = ^ + m( R ^L + ^\ iflS, 



c 2 V dR dz I c 2 dR c 2 dz 



(84) 



where the dot denotes derivation with respect to t. In particular, equatorial circular orbits 
must satisfy the conditions R — z — 0, R — z — and z — 0. In this case ( 152"]) reduces to 



which can be used to derive an expression for v 2 at 1PN order. In particular, note that we 



can use the Newtonian expression t> 2 = Rd<p/dR each time that v 9 is accompanied by an 
inverse power of c. After some straightforward algebra we get 



where it is understood that all derivatives are evaluated at z — 0. There is a crucial difference 
between the above relation and the classical formula for the rotation curves: in the Newtonian 
case is linear in d<j)/ dR, whereas in the 1PN case, it depends in nonlinear terms involving 
and derivatives of the gravitational potentials. This nonlinear dependence may be significant 
in some cases. 
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B Virial Theorem in the 1PN Approximation 

The virial theorem is an important general result of the kinetic theory relating the contribu- 
tion of kinetic and potential energy to the temporal change of the inertia tensor associated 
to the system. In order to obtain the 1PN version of the virial theorem we can start from 
the conservation laws 



dT ^u 
dx^ 



i7) 



At first order in v/c, this expression reproduces the Newtonian mass and momentum con- 
servation laws, which are given by 



1 d T 00 d T i0 
c dt dx i 

1 2 
ldT 0i Q rpij 

+ 



c dt 



dxi 



1 dt ^oo 



c 2 dx % 



(88) 
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The corresponding 1PN corrections of the above laws are obtained by taking into account 
Christoffel symbols at different orders [5T] 



2 3 

1 d T 00 d T i0 

c dt dx 1 

3 4 

1 d T 0i d T ij 
+ 



c dt 



dxi 



c 3 dt 

- J&T 00 - J&T 00 -(2 kj ko kk) T°> 
-(k +5 lk h +5 ik k) k 



d (2(j) 2 iff 



.dx 1 V c 4 



I%1 > 
c 4 dt 



c 2 dx' 



2 

rpOO 



(90) 



dx* 13 dt) A^dx* dlk dxi) A } 



On the other hand, the energy-momentum tensor, given in ( fl71) . can be split in two first-order 
contributions: 



rpflU 



F V^V u d 3 v + 



o /3v 2 
F 



2 

+ F 



V»V u d 3 v 



(92) 



where = (c, v). After some definitions 



o 



/" , 2 /■ 2 „ 2 2 / 

/ F d 6 v, P= F d v, and p = p + pl 



o /3v 2 



(93) 
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and using the different probability densities to compute expectation values, 

— o i f c 



2_ 

A 



or 1 



o /3v 2 60 



Adrv. 



(94) 



we get the following simplified expressions for the momentum-energy components at different 
orders: 



(95) 
(96) 






2 


2 

y00_ 


c 2 p, 


1. 0- 

T 0l = cpv\ 


3 

rpOi 


2^ 

c Pv l , 


2 

rpij_ 




:/5 V l Vl , 


4 2 2 

T 3 =pv^ 



Then, equations fl90l) - fl9Tl) . can be written as 

d(|) g(|g) P_30 
dt dx i c 2 dt 



(97) 



2 A 



d(pv l ) d{pv i vi) 



dt 



dxi 



- P 



dx l V c 2 c 2 / c 2 <9i - 



2 

- P 



c 2 V&r- 7 9a: 4 



c 2 V jk dx* 



ik 



dxi 



).(9* 



Now we can compute the symmetrized integral of the product x k /2 times the above equation: 



where 



2 



K 



ik 



1 d 
2dt 



2 , J_ 



2 2 



cf a; p v i +x l v k ) = 2 K lk + W 



(99) 



— / d 3 xf> x' 



c 2 <9t c 2 <9x J 



ft 



+ 



c 



2 A ,^-u„^ 



' dxi ) 



-- I d 3 x P x' 



dx k \c 2 c 2 / c 2 dt c 2 \dxi dx k ^ dt 



+ 



c- 



45, 



ki 



dxi J 



(100) 
(101) 
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By defining 

2 



2 



I ik = J d 3 xpx l x k , (102) 
we can write the time derivative of the equation ( 190]) as 

2 ik ~ 2 2 

i!L = J d 3 x^£x i x k = J d 3 x 2 p (x k v* +x l v*) + j d 3 x^^x l x k (103) 
and its second time derivative 

d 2 P k d f , s 2 / 2 ^ rf , 



^ [ d 3 xp {x k v* W v k ) + 4- [ d 3 x^ d 4-x l x k . (104) 
dt J H y 1 dt J c 2 dt y ' 



dt 2 dt 

Putting all together we obtain the tensor post-Newtonian virial theorem 



2 o 



J2 Tik 2 2 1 d f P 

2K ik + W ik +-— d 3 x^^x l x k , (105) 



dt 2 2 ^ 7 c 2 dt 

and taking the trace, we obtain the scalar post-Newtonian virial theorem 

d 2 I 2 2 1 d f p 90 9 , , 

where 

2 2 -, 

0= trace(O lfc ) . (107) 

In summary, we can collect all the above results and state that if we define a moment of 
inertia tensor, 

lk - I I ' 1 + ^ - ^ ] x l x k Fd 3 vd 3 x, (108) 




a kinetic energy tensor, 




3y2 6< M:zr^3 ,3 



K lk = - / 1 + — £ v i v k Fd i vd 6 x, (109) 



2./ ./ I" ' c 2 c 2 
and a potential energy tensor 
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where is the Newtonian potential energy tensor and W is given by fllOip . then they 
satisfy the following relation 



o 



<]1 '''' 2K lk + W ik + ~ [ d 3 x^x l x k , (111) 



dt 2 2 dt J c 2 dt 

which can be enunciated as the 1PN tensor virial theorem. The scalar virial equation is 
obtained by taking the trace of the above relation: 

rf2/ nr^ y t t- 1 d f j3 P 8<P 2 

^ = 2K + W+ 2dt 



C Jacobi's Integral in the 1PN Approximation 

For the purposes of the present paper it is important to have in hand an expression of the 
Jacobi's integral at 1PN order, in order to verify that it has the same form as in Newtonian 
theory, i.e. Ej = E — QL Z . 

In a rotating reference frame with angular velocity fl, velocities are related through 



v — fl x x, 



(112) 



where v and v are the velocity measured from the rotating and inertial frame, respectively. 
From this relation, and the 1PN equations of motion it is possible to derive a 1PN corrected 
version for Ej [65J: 



Ej 



d 2 



1 <~ ,2 I 



+ (£1 x x) 2 + 



3v 4 2 , 3<f)v 2 , n . . 30. ^ , 2 
J + y + ^y-(^x)^ + y(Oxx) 2 



-i(n x x) 4 + j(n x x) 2 + v • (n x x )5 2 + ~(v • (n x x)) 2 

In particular, if O = i?e z and x = Ren, the above formula reduces to 



• (113) 



Ej = - + 



Q 2 R 2 _L_ 
2 + 2c 2 



r.l 



3v 

— + 2 + 2^ - 3# 2 - 2QB£ V + 30i7 2 i? 2 



i7 4 _R 4 v 2 f2 2 R 2 „ _ _o „o 

— — + + 2QRv^v 2 + v%Q 2 R 2 



Now, by replacing v R = Vr, v z = v z , and v v = v v — QR in (I114p . we can write 

Ej = E — f2L z , 



;ii4) 



(115) 



where E and L z are given by f fTUj) and ffTTj) . respectively. Thus we see that the total energy 
associated to a test particle, measured from a rotating frame, has the same form as in the 
Newtonian theory. 
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D Derivation of the Constants Bm and Cm 

For the DF given in ( 167|) and for any function S{R 1 Vr 1 v v ) 1 we can write 



Sfd 2 v = 2^ 4 k / dv v / 



(Z,2 _ „2 )l/4 



(116) 



where 



45vrGM 



/i = V^ei^ 1 ^' ^ = V^ 2 -^-^) 2 - ( 117 ) 

In particular, by setting S = v v in (I116P and introducing the result in fl60|) we obtain 



8 V 2vra 5 

n=0 



p r ^ 5 / 15GM3 /i 2 4 



from which can be found by using the orthogonalization condition J^_ 1 P m (x)Pi(x)dx 
25 lm /(2l + l). We find 



5 hr>(?M 3 r 1 , 

C 2n = (4n + l)|y-^g- jf v / T^x 4 P 2n (x)dx. (118) 

The above integral can be solved by setting x = cos 9 and using the following identity: 
tt/2 TrrrfHi+l)! 2 

(sin^) m P 2n (cOS^)^ = ; r 2 { , r— Ti r. (119) 

o 1 ' K ' 2r 1 + n + f r | + f -n r 1+n r |-n 1 ; 



After some calculations, we get 



tt 2 (45 - 48n - 92n 2 + 16n 3 + 16n 4 



Vl - x 2 x A P 2n (x)dx = K —, ^Tv r-^ r— ; (120) 

v 7 i28r(| - n)r(| - n)r(i + n)r(4 + n) 

Plugging the above result in (11181) we obtain the relation (!77j) . 

Now, performing the same procedure on equation (1591) and obtain 



2^ r, / \ 757r 2 / GM\ . 32t/> 160irGM , 



n=0 

where 



H{^) = 375r7f - 154r/, 6 + 392r/ 4 + 14r/ 2 - 217. (121) 



Again, we use the orthogonality of Legendre polynomials to provide an explicit relation for 
B 2n . In the calculations we also find useful (i) the formula 

1 2mD / \ t v^2- 2m T(l + 2m) 

x 2m P 2n (x)dx 



-i 



r(l + m — n)r(^ + m + n) 
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and (ii) the expression 

/l m m—k 

Vl - x 2 P 2n (x)P 2m {x)dx = ^2^2 
1 z — n A— n 



7r(4m - 2k)\(m - k)\ 



2 2m k\{2m - k)\j\(m - k - j)!(2m - 2k)\ 

(zM^±ia! (122) 

r(2+i + n)r(j-n + |)r(i + n)r(i-n)' v ; 

which we derived from (11191) and the series representation of the Legendre polynomials. After 
some calculations we get 



y5 2n92n (0) + 2n 2 aGB 2n = (^f) 



2 



oo 

^2n ~ n2n,2m 



m=0 



where 



with 



75tt 5 / 2 (4n + l)J2l =0 a k n k 
2n 16384 r(5 - n)r(^ + n) 



(123) 



(124) 



and 



a = 370224 ai = -337854 a 2 = -628841 
a 3 = 185300 a 4 = 174491 a 5 = -25768 

a 6 = 16600 a 7 = 992 a 8 = 496 (125) 



257T 3 Z" 1 

n 2 n,2m = J Vl - X 2 P 2n (x)P 2m (x)dx. (126) 



From (I123p . it is straightforward to derive ( !7H|) . 
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